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The time-dependent superfluid local density approximation (tdslda) is an extension of the 
Hohenberg-Kohn density functional theory (dft) to time-dependent phenomena in superfluid 
fermionic systems. Unlike linear response theory, which is only valid for weak external fields, 
the tdslda approach allows one to study non-linear excitations in fermionic superfluids, including 
large amplitude collective modes, and the response to strong external probes. Even in the case of 
weak external fields, the tdslda approach is technically easier to implement. We will illustrate the 
implementation of the tdslda for the unitary Fermi gas, where dimensional arguments and Galilean 
invariance simplify the form of the functional, and ab initio input from quantum Monte Carlo (qmc) 
simulations fix the coefficients to quite high precision. 



Iinear response theory is a popular tool for study- 
ing the dynamics of a quantum many-body system. 
_ Formally, the change in the number density (often 
referred to as the transition density) in response to a 
weak external potential V ext (f,t) is given by 



6n(f,t) 



df'dt'n(f,t,f',t')V ext (f',t') 



(1) 



where TT(r, t, f',t') is the linear response function of the 
system. Since, for a system in equilibrium, U{f, t, f', t') 
depends only on the difference t — t', one usually works 
with the Fourier transforms: 



6n(f, a>) 



df'n(f,f',u>)Vext(f',a>) 



(2) 



The linear response function TT(f, f , tu) has poles at fre- 
quencies corresponding to the various excited states of 
the system, which allows one to express these excited 
states in a form independent of the external probe: 



df'A(f,f',a>)5n(f',tu) = 0. 



(3) 



Here A{f,f',w) is the operator inverse of TT(f,f', tu). 
The existence of A(f, ?', cu) is nontrivial as the operator 
TT(r, f , a>) may be singular due to zero modes (Goldstone 
modes) arising from various conservation laws. 

This approach is appealing, because solutions of 
Eq. describe intrinsic excitations of the system. How- 
ever, it is clearly limited to describing small amplitude 
excitations where the response remains linear and the 
external potential is weak. These equations are also 
technically difficult to solve due to the high dimension- 
ality of the matrices involved: especially in the case 
of inhomogeneous systems. This makes it practically 
impossible to implement a fully three-dimensional cal- 
culation, and they have only been solved in systems 
with a high degree of symmetry: infinite homogeneous 
systems for example, or axially / spherically symmetric 
configurations. Even in such cases, limiting assumptions 
or approximations are often required. 



Here we shall describe a different approach: time- 
dependent density functional theory (tddft). This not 
only allows one to study non-linear excitations, but also 
allows one to consider fully three-dimensional equations. 
Although exact in principle, there is no simple prescrip- 
tion for computing an exact density functional in a non- 
perturbative theory (see [1] for recent discussions), and 
one must first formulate an approximate functional that 
captures the relevant physics. In the case of the unitary 
Fermi gas, the lack of scales greatly restricts the possible 
forms for the functional, and an extremely simple form 
— the superfluid local density approximation (slda) [2] 
(described in Sec. |I A) — appears to capture much of the 
relevant physics. The time-dependent superfluid local 
density approximation (tdslda) requires one to solve 
a system of coupled time-dependent three-dimensional 
nonlinear Schrodinger-like equations of the form 



ih 



at 



[H(f,t)+V ext (f,t)]¥ k (f,t) 



(4) 



Here W^(f,t) is a vector of single-quasiparticle wave- 
functions, the exact meaning of which will be explained 
below, and the corresponding single-particle Hamilto- 
nian H(f, t) is a partial differential operator. The main 
complexity of this system of equations arises from the 
fact that the single-particle Hamiltonian H(f, t) depends 
non-linearly on all the single-quasiparticle wavefunc- 
tions W^{f,t). The simplification is that H contains only 
differential operators (no integral operators either in 
time or space), and can be efficiently applied on each 
wavefunction independently, allowing the method to be 
efficiently parallelized. Since no matrix operations are 
involved (the kinetic and potential parts are applied sep- 
arately and efficiently using the fast Fourier transform 
(fft), and the memory requirements are significantly 
reduced compared to solving Eq. ([3 1. 

The tdslda also has conceptual advantages over some 
traditional approaches to superfluid dynamics: unlike 
two-fluid hydrodynamics, the tdslda can correctly de- 
scribe quantized vortices and their dynamics, and con- 



2 



tains naturally the critical flow velocity at which a super- 
fluid can turn into a normal fluid; in contradistinction 
to the Gross-Pitaevskii or Ginzburg-Landau approaches, 
the normal fluid to superfluid transition is within the 
scope of the theory. Moreover, a number of large am- 
plitude collective modes have been studied with the 
tdslda that defy a description within two-fluid hydro- 
dynamics, Ginzburg-Landau, or Gross-Pitaevskii frame- 
works 



in Sec. I A The technical challenges are: 1) diagonaliz- 
ing the single-particle Hamilton (|6j; 2) solving the self- 
consistency equations to determine stationary (ground 
state) configurations; and 3) stably and self-consistently 
evolving the single-particle states Q to describe the dy- 
namics. Typically one applies all three techniques, first 
solving for an initial stationary configuration, then driv- 
ing the system to explore the dynamics — stirring to 
generate vortices for example. 



I. METHODOLOGY 



A. The Functional 



A precise formal statement of a density functional the- 
ory (dft) starts with some physically motivated en- 
ergy functional E[ni , n.2, ■ ■ • ] of various densities n^(f, t). 
To simplify the formal structure, we express this as a 
function of the density matrix E(p) though in the end 
we shall only consider local functions (see Sec. |I A) . 
Once specified, one simply minimizes the free energy 
F(p) = E(p) +TTr(plnp) subject to the normalization 
constraint on p + Cp T C = 1 dictated by Fermi statistics, 
where C = C T is the charge conjugation matrix. The 
constrained minimization of the functional F(p) results 
in the standard Fermi distributiorQ 



p = 2jk)n FD (E k )(k| 



1 +e P(H(p)-CHT(p)C)' 



to obtain the following equations of motion 

H(p)|k) = 6 ^ ) |k)=E k |k>, 
5p 

p = 2jk)n FD (E k )(k|, 



(5) 



(6a) 
(6b) 



which must be solved self-consistently. The eigenvalues 
E k are the Lagrange multipliers of the associated normal- 
ization constraint. The formulation of the tddft follows 
simply by using H(p) to generate the time evolution of 
the single particle states, 



ih9 t |k) =H t (p)|k), 



(7) 



typically in the presence of some time-dependent ex- 
ternal potential included in H t (p) = H(p) + V ex t(t), for 
example, or a gauge coupling in the case of an electro- 
magnetic external field. 

The physical content of the dft enters through the 
formulation of the function E(p) as we shall discuss 



In practice, one does not work explicitly with the density 
matrix p but rather with a set of physically motivated 
local densities. It is convenient to express these concepts 
in the language of second quantization. We consider 
two species with operators c-|- and cj_ representing two 
hyperfine states. 

For a two component system, the most general wave- 
function that allows for all possible pairings has four 
components: W = (c^,C|_,cj,c|). In terms of components 
of the wavefunction, we will write HW^ = E^^ where: 
^k(r,t) = (f|k) = (u kt (f,t),u u (f,t),v kt (f,t),v u (f,t)). 
In what follows we shall drop the explicit (f, t) depen- 
dence. In this formulation, the time evolution of a single- 
particle wavefunction ¥ k is: 



3t 



v kt 

/h t + U t 





A* 



X 

H^ + U 
—A* 




■LL 



A 


-x* 



Uk4 
Vkt 



(8) 



where h^j, = —V 1 /{2m^^), 11-^ is the self-energy, and 
A oc (cjcj) is the pairing field. One needs this full 
four-component formalism if x <x (cjcj_) 7^ 0. (A spin- 
orbit coupling in the nuclear problem would appear 
here for example.) For the unitary gas, however, we 
consider only attractive s-wave interactions (thus, x = 0), 
allowing us to express everything in terms of the usual 
two-component Bogoliubov-de Gennes (BdG) form ¥ k = 
(u k ,v k ): 



u k 

"at V v k 



h-j- + LU 
A* 



U 



u k 
v k 



(9) 



1 Formally, this constraint can be implemented using a Lagrange 
multiplier, but it is much easier to see the results by letting p = 
1 /2 + x — Cx T C where x is unconstrained, and then performing 
the variation with respect to x. 



Note that the structure of these equations is that of a 
single quasiparticle Hamiltonian: Indeed, for the choice 
of functional we consider below, this will look formally 
like the standard Bdc equations, however, the coeffi- 
cients will be determined from the functional rather 



3 



than from a direct mean-field approximation of a mi- 
croscopic theory. In the usual formulation of a dft for 
normal systems, the single particle states need not bear 
any formal relationship to the physical quasiparticles. 
Within the slda, however, we have found that the quasi- 
particle properties — their dispersion relationship for 
example — can also be successfully modeled with the 
appropriate choice of functional. 

For simplicity we shall consider here only the sym- 
metric case tvi- = Tt|^ = n + /2 where the two states have 
identical masses and describe the slda. (See [4! for de- 
tails about the asymmetric slda (aslda) extension.) We 
consider three densities and one current: 

n+(f) =2^iv k (f)| 2 n FD (-E k ) ~ (10) 
x+(f) =2^[Vv k (f)| 2 n FD (-Ek) ~2jVc+ a -Ve ff >, 
v(f) = \Y_ u v (f)v k (f) [n FD (-E k ) - n FD (E k )] - (c t c 4 ), 

k 

) + {f)=iY_ [v k (f)Vv k (f)-v k (f)Vv k (r)] n FD (-E k ). 

k 

We use the kinetic energy density t+ in the spirit of 
Kohn-Sham, and the anomalous density v to account for 
pairing within a local theory. For time-reversal invariant 
ground states, the current density j + vanishes. It must 
be considered when considering time-dependence to 
ensure that the energy density is covariant under local 
Galilean transformations. In nuclear physics Galilean in- 
variance have been considered for quite some time (flEl, 
and the contribution of these currents is often essential 
for describing the properties of excited states. It is easily 
demonstrated (see [4] for details) that when changing to 
a frame with velocity v, the currents and kinetic densities 
transform as 

j*+ -> f+ + Mvn + , t+ -> t+ +v ■ j* + + 2M|v| 2 n + (11) 

where M = JVLf = M± is the bare mass of the particles. It 
follows that for symmetric two-component systems, the 
following is Galilean invariant: 



IT + ! 2 

2Mn, 



(12) 



Fhe center of mass motion may be separated from the 
intrinsic energy density (the total energy E = Jd^3 ]f£) 



D+r 

2n+ 



+ £(T + ,n+,v) 



(13) 



and the two always enter the functional as 



n\ /3 /y-A/<x.' 



(14) 



where y parametrizes the pairing strength, <x = M./M e g 
is the inverse effective mass, and A is a momentum space 
cutoff. The most straight-forward functional constructed 
from these quantities is the slda: 



£slda(t + ,ti + ,v) 



h 2 

M 



y'y 



n! /3 /T-A/a 



+ 



+ pA( 37t 2 ) 2/3 n 5/3 ) ^ (15) 



Varying this functional leads to the following identifi- 
cation of the single particle Hamiltonian h = ru = h|_, 
potential U, and gap parameter A: 



H : 
A : 



-0L- 



2M 



n\ /3 /y-A/a' 



2\2/3 2/3 

1 ' n / 



At A 

3T^ 



+ V, 



ext- 



(16) 
(17) 

(18) 



For spatially varying systems, momentum is not a good 
quantum number and a simple momentum space cutoff 
cannot be implemented. Instead, one can use an energy 
cutoff, limiting the sums in Eqs. ( [io| for energies |E k | < 
E c . The homogeneous equations can then be used to 
translate this into a position dependent A(f) that may 
be used in the previous equations and which has very 
good convergence properties [9] (see also M): 



Mk c (f) 



ko(f) . k c (f)+ko(f)' 

2Mf)\m-kom.^ (19) 



where ko and k c are defined by 



a- 



h 2 k 2 (i 
2M 



M.+ U(f)=0, oc 



2M 



^ + U(f) = E c . 



To complete the functional, one must determine the 
parameters oc, |3, and y. We do this by matching the pre- 
dictions of the functional in the thermodynamic limit to 
accurate quantum Monte Carlo (qmc) calculations. Fit- 
ting the energy and quasiparticle spectrum determines 
the following values for the unitary gas (see [4I for a 
detailed discussion of this fitting procedure): 



such that £ is locally Galilean invariant. 

The form of the functional is further restricted by the 
fact that the anomalous density v(f,f) ~ (c-|-(f)cj,(f')) oc 
|f — f'| _1 is ultraviolet divergent in the local approxima- 
tion. This divergence also appears in the kinetic term t + 



a =1.094(17), [3 = -0.526(1 i 



-0.0907(77). 



The tdslda satisfies all expected conservation laws: 
energy in the absence of time-dependent fields, lin- 
ear/angular momentum if the corresponding symme- 
tries are not broken, gauge and Galilean invariance, and 
particle number in the absence of applied external pair- 
ing fields. 
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B. Technical Notes 

Solving the self-consistency conditions requires solv- 
ing such a large number of simultaneous equations that 
typical root finding methods employing a Jacobian com- 
putation are prohibitive. However, treated as an iterative 
method — take an initial set of densities, form the po- 
tential \l6) , diagonalizing the Hamiltonian to obtain a 
new set of single particle wave functions, and then form 
a new set of densities \io\ — the self-consistency cycle 
is typically close to convergent. As a result, a mem- 
ory limited implementation of Broyden's method ffToll 
works well to accelerate convergence, thereby determin- 
ing equilibrium configurations to use as an initial state 
for a subsequent time-dependent simulation. 

The output of this is a complete set of wavefunctions, 
typically represented on a periodic lattice. These are 
then fed into the time dependence equations Q to gen- 
erate the time-dependent states |n(t)). Note that at each 
time-step, the Hamiltonian must be updated to reflect 
the current ensemble of states. We have found that a 
multistep predictor-modifier-corrector method due to 
Adams-Bashforth-Milne (see [11]) works well (see HT2I1 
for implementation details and parallel scaling perfor- 
mance.) Periodic lattices enable us to use the fft to 
efficiently transform the wave functions between po- 
sition and momentum space so that the kinetic and 
potential parts of the Hamiltonian may be applied by 
simple multiplication. This allows us to perform fully 
three-dimensional simulations. 

C. Validity domain 

If one has an exact density functional, then the tddft 
technique can be shown to deliver the exact time- 
evolution of the one-body density Ji3] - [T5| . If one is 
interested in higher-order operators, however, then ex- 
tensions to the technique are required |i6J. These are 
significantly more costly, but still within computational 
reach for carefully chosen problems. 

The main limitation is that an exact density func- 
tional is not known. Thus, the dft requires careful 
benchmarking to determine the domain of validity. At 
present, the slda has been formulated and fit using 
QMC calculations of the T = thermodynamic limit of 
the three-dimensional unitary Fermi gas. This has been 
benchmarked against trapped systems to an accuracy of 
a few percent [4], indicating that the omitted gradient 
corrections are quite small. Thus, the slda is reliable 
for cold symmetric systems up to small gradients cor- 
rections. The asymmetric extension (the aslda) has also 
been formulated and fit to QMC data. The extension to 
finite-temperatures is still an open problem. 



D. Relevance to other theories 

The aslda subsumes the usual mean-field Bdc equa- 
tions, but extends these considerably. For example, it 
includes a self -energy contribution that is neglected in 
the zero-range limit of the mean-field Bdc equations. 
The aslda lacks the variational property of the mean- 
field BdG equations, but with careful validation, has the 
ability to provide a much more quantitatively accurate 
description of fermionic superfluids [gj. 

II. APPLICATIONS 

We present here briefly two quite spectacular results 
obtained using the tdslda in a unitary gas. We prepare 
a system in its ground state in an axially symmetric trap 
(with an essentially flat bottom) and homogeneous with 
periodic boundary conditions in the third direction. We 
then adiabatically introduce two types of quantum stir- 
rers: 1) a spherical projectile flying along the symmetry 
axis with a speed v p = 0.2 vp (where vp is the Fermi 
velocity); and 2) a rod parallel to the symmetry axis with 
a diametrically opposed sphere (breaking translational 
invariance along the tube) moving with a constant an- 
gular velocity about the center of the tube and a linear 
velocity lower than the critical velocity of the unitary 
Fermi gas v c w 0.365 vp J17) [i8l . In the first case, the 
spherical projectile, after passing through the system, 
generates a rather elusive excitation mode of a super- 
fluid: a vortex ring. In the second case, the two quantum 
stirrers (the rod and the sphere) generate five vortices. 
The sphere breaks the translational symmetry, exciting 
Kelvin modes along the vortices, and, at the same time, 
exciting phonons in the superfluid to form a compli- 
cated mixture of dynamical modes. In each of these 
simulations we solved about 22000 time-dependent 3D 
coupled nonlinear partial differential equations on a 32 3 
spatial lattice for a sufficiently long period of time. 

III. RELEVANCE TO OTHER SYSTEMS 

Even though we have only illustrated the tdslda in 
the case of a unitary Fermi gas, this is a rather general 
approach suitable to describe the dynamics of virtu- 
ally any fermionic superfluid with s-wave pairing. The 
tdslda has already been used to describe nuclear sys- 
tems: in particular, the first attempt to describe induced 
nuclear fission was recently performed. Although not 
yet explored, it appears that the extension to pairing in 
other partial waves (p-wave and d-wave for example) is 
straightforward. 



Figure 1. Two frames of 3D time dependent simulations of a unitary Fermi gas confined to a cylindrical trap and subject to a time 
dependent external potential. On the left, a hard sphere moved along the trap axis, generating a vortex ring in its wake. On the 
right, the external potential was a vertical rod and a diametrically opposed sphere which stirred the system, generating five 
vortices. Kelvin waves have been excited along each vortex. The last two vortices have been generated simultaneously: they are 
essentially on top of each other and separate at a later time. 
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